//
// oversampling a/d
//
// for a 16 MHz clock, the adc rate is 9615.4 sps and the 512-stack is
// 18.78 sps.  the final sample rate is our only assurance that we
// haven't skipped interrupts.
//
// we actually use a 1024-long window stepped at 512.  using the
// double window should (i think) put the rectangular window's first
// spectral zero at the nyquist frequency for the final sample rate.
//
// 24 July 2011  went to 4 x 512 windows.  the aggregate window is 2048 points
//   and is stepped at 512.  the first zero should be at 18.78/4 or about
//   4.7 Hz.
//
// worth noting that we use the 512 window to set the output sample rate
// while the number of 512s that are combined scales the spectral response.
//
// 24 September 2011  added digital filtering.  also tried to minimize
//   start-up transients.

#include <avr/io.h>
#include <avr/interrupt.h>

volatile unsigned long runningsum;
volatile unsigned short runningcount;
volatile unsigned long prev3;
volatile unsigned long prev2;
volatile unsigned long prev1;
volatile unsigned long current_sum;

volatile boolean next_sample_ready;
volatile unsigned int next_sample;

//
// Uncomment one of the under/define pair below.  If FILTER is
// undefined the program implements a benign length-one allpass
// filter.  If FILTER is defined the program (currently) implements a
// 55-coefficient lowpass filter with a cutoff of 2.5 Hz.
//
// Because of a deficiency in the arduino ide we have to copy the
// filter coefficients directly into the code below.
//
//
#define FILTER

#ifdef FILTER
//
// If you want to replace the righteous filter below, just replace the
// coefficient list.
//
// ScopeFIR design: simple lowpass 1.8-2.6, stop 60,  59 terms
//
const double coeff[] = {
  0.003288038758116306,
  0.011554564141729728,
  0.028285173151443389,
  0.055293154901223225,
  0.091740492086433742,
  0.132882755087075030,
  0.170180567155164610,
  0.193144482856034080,
  0.192529777511811050,
  0.163937396952292090,
  0.110463646637543680,
  0.042818962194838196,
  -0.023110756249503973,
  -0.071124795415197492,
  -0.090111032202806640,
  -0.077964180145652706,
  -0.042326982981975958,
  0.002170770746045263,
  0.039266660956168702,
  0.056858755659304719,
  0.051139515596325083,
  0.027313694379628065,
  -0.003195962263591410,
  -0.027698476190158871,
  -0.037344306211302712,
  -0.029992980672969630,
  -0.010491584708780343,
  0.012199628246658469,
  0.029051229108537493,
  0.034558716072542808,
  0.028404241341067425,
  0.014750702617246883,
  -0.000087250486180852,
  -0.010558388571706344,
  -0.013936832803067312,
  -0.010863198142809184,
  -0.004330600742312823,
  0.002067539687054484,
  0.005769049959253918,
  0.005970771233240058,
  0.003625233538429873,
  0.000480788679567260,
  -0.001876370459966127,
  -0.002672068267129423,
  -0.002072029155244571,
  -0.000819126898662612,
  0.000384410977839965,
  0.000825152340957656,
  0.001180138997267963,
  0.000378073825107134,
  0.000034272719945145
 };
#else
//
// default allpass filter
//
const double coeff[] = {
  1.00000e+00,
};
#endif

const unsigned int ncoeff = sizeof(coeff) / sizeof(coeff[0]);
unsigned int lagarray[ncoeff];
unsigned short first_adc = 10; // initial samples to skip
unsigned short first_loop = 1; // deal with the first composite sample


void initADC() {

  // don't depend on global initializers
  runningsum = 0;
  runningcount = 0;
  prev1 = prev2 = prev3 = 0;
  current_sum = 0;

  next_sample_ready = false;
  next_sample = 0;

  // internal AVcc ref, no left adj, channel 0
  ADMUX = _BV(REFS0);
    
  // ACME off, free-running mode
  ADCSRB = 0;

  // enable ad, intr, auto
  // set divisor to 128
  // start first conversion
  ADCSRA = _BV(ADEN) | _BV(ADIE) | _BV(ADATE)
    | _BV(ADPS0) | _BV(ADPS1) | _BV(ADPS2)
    | _BV(ADSC);
}


ISR(ADC_vect) {
  byte low = ADCL;
  byte high = ADCH;
  if(first_adc) {
    --first_adc;
    return;
  }
  runningsum += low + (high << 8);
  if(++runningcount == 512) {
    if(first_loop)
      prev3 = prev2 = prev1 = runningsum;
    else {
      prev3 = prev2;
      prev2 = prev1;
      prev1 = current_sum;
    }
    current_sum = runningsum;
    // if we sum two we have to unshift 4
    // if we sum four we have to unshift 5
    next_sample =
      (unsigned int) ((current_sum + prev1 + prev2 + prev3) >> 5);
    runningsum = 0;
    runningcount = 0;
    next_sample_ready = true;
  }
}


void setup() {
  unsigned int i;
  Serial.begin(9600);
  Serial.flush();
  first_loop = 1;
  first_adc = 10;
  initADC();
}


void loop() {
  double fsum;
  unsigned int i;
  unsigned int filtered_signal;

  if(next_sample_ready) {
    next_sample_ready = false;

    if(first_loop) {
      // flood the array with the first value
      for(i = 1; i < ncoeff; i++)
	lagarray[i] = next_sample;
    } else {
      // put the new sample into the bucket brigade
      for(i = ncoeff - 1; i > 0; i--)
	lagarray[i] = lagarray[i - 1];
    }
    lagarray[0] = next_sample;
    first_loop = 0;

    fsum = 0.0;
    for(i = 0; i < ncoeff; i++)
      fsum += lagarray[i] * coeff[i];

    filtered_signal = (unsigned int) (fsum);

    Serial.println(filtered_signal, DEC);
    Serial.flush();

  }
}
